%plots perturbation population response
%plots population resonse on running onset (fb+dark)

ii=figure;
oo=figure;
perturb_act=[];
run_act=[];sem_act=[];
l_win = 10;
r_win = 50;
th_run = 0.01;
crr=ones(6,10)*NaN;
colors=[0 0 1;0 0 1;1 0 0;1 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;0 0 0;];
plotit=1;
for site=1:size(proj_meta,2)
no_time_points = size(proj_meta(site).rd,2);

act=[];
type_sessions={'f','p','d'};
for tp=1:no_time_points
    act=[];
    for lyr=1:size(proj_meta(site).rd,1)
        act=[act; proj_meta(site).rd(lyr,tp).act];
    end
    
    sess_start=proj_meta(site).rd(lyr,tp).nbr_frames;
    sess_start=cumsum(sess_start);
    sess_start=[1 sess_start+1];    
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(1))==1);
    indices=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    
    velM = proj_meta(site).rd(lyr,tp).velM_smoothed;
    velP = proj_meta(site).rd(lyr,tp).velP_smoothed;
    ps=proj_meta(site).rd(lyr,tp).ps_id;
    run=velM/max(velM)-mean(velM);
    
    velM_r = run>th_run;
    start_perturb=find(diff(ps)>1);
    start_perturb=intersect(indices,start_perturb);
    start_perturb(start_perturb<=l_win)=[];
    start_perturb(start_perturb>=length(ps)-r_win)=[];
    ind_perturb=[];
    for kkk=1:length(start_perturb)
        if sum(velM_r(start_perturb(kkk)-l_win:start_perturb(kkk)+r_win))>ceil((l_win+r_win)*.9)
            ind_perturb = [ind_perturb start_perturb(kkk)-l_win:start_perturb(kkk)+r_win];
        end
    end
    
    %     figure,temp3=mean(act(:,indices))/max(act(:));plot(temp3),hold on, temp4=velM;temp4(setdiff(1:length(temp4),ind_perturb))=0;temp4=temp4(indices);temp4=temp4+min(temp3);plot(temp4,'r')
    temp2=act(:,ind_perturb);
    temp2= reshape(temp2',l_win+r_win+1,numel(temp2)/(l_win+r_win+1))';
    perturb_act(site,tp,1:l_win+r_win+1)=mean(temp2,1);

    %% running onset
    % first fb session - filter out the running and stim moving
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(2))==1);
    indices=[];
    ind_run=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    velP_r = velP/max(velP)-mean(velP);
    velP_r = velP_r<th_run;
    temp_start_run = velM_r(1:end-2)+velM_r(2:end-1)+velM_r(3:end);
    temp_start_run(temp_start_run>1.5)=3;
    temp_start_run(temp_start_run<=1.5)=0;
    temp_start_run = strfind(logical(temp_start_run),[0,1]);
    start_run=intersect(indices,temp_start_run);
    start_run=intersect(start_run,find(velP_r==1));

    for kkk=1:length(start_run)
        if start_run(kkk)>10 && start_run(kkk)+30<sess_start(end)-1
        if sum(velM_r(start_run(kkk):start_run(kkk)+30))>28 && ...
            sum(velP_r(start_run(kkk):start_run(kkk)+30))<28
                ind_run=[ind_run start_run(kkk)-l_win:start_run(kkk)+r_win];
        end
        end
    end
    act2=mean(act(:,indices));
    act2=act2/max(act2)-mean(act2);
    crr(site,tp)=diag(corrcoef(act2,run(indices)),-1);
    % for dark
    ac_sess=find(strcmp(proj_meta(site).rd(lyr,tp).session,type_sessions(3))==1);
    indices=[];
    for kk=1:length(ac_sess)
        indices=[indices sess_start(ac_sess(kk)):sess_start(ac_sess(kk)+1)-1];
    end
    start_run=intersect(temp_start_run,indices);

    for kkk=1:length(start_run)
        if start_run(kkk)>10 && start_run(kkk)+50<sess_start(end)-1
        if sum(velM_r(start_run(kkk):start_run(kkk)+30))>28
            ind_run=[ind_run start_run(kkk)-l_win:start_run(kkk)+r_win];
        end
        end
    end
    temp2=act(:,ind_run);
    temp2= reshape(temp2',l_win+r_win+1,numel(temp2)/(l_win+r_win+1))';
    run_act(site,tp,1:l_win+r_win+1)=mean(temp2,1);
    

%%
perturb_act(site,tp,:)=perturb_act(site,tp,:)/mean(squeeze(perturb_act(site,tp,1:10)));
run_act(site,tp,:)=run_act(site,tp,:)/mean(squeeze(run_act(site,tp,1:5)));
if plotit
    figure(ii);
    % set(0,'DefaultAxesColorOrder',colors(1:no_time_points,:));
    subplot(ceil(size(proj_meta,2)/2),2,site);
    hold on
    plot(squeeze(perturb_act(site,tp,:))','Color',colors(tp,:))
    % plot(perturb_act');
    title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
    
    figure(oo);
    % set(0,'DefaultAxesColorOrder',colors (1:no_time_points,:));
    subplot(ceil(size(proj_meta,2)/2),2,site);
    hold on
    plot(squeeze(run_act(site,tp,:))','Color',colors(tp,:))
    % plot(perturb_act');
    title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end
end
end

% average tp - mismatch
end_points = [5, 6, 9, 9, 9, 10];
figure;
for site = 1:size(proj_meta,2)
    subplot(ceil(size(proj_meta,2)/2),2,site);
    hold on
    plot(mean(squeeze(perturb_act(site,1:2,:)))','b')
    plot(mean(squeeze(perturb_act(site,3:4,:)))','r')
    if end_points(site) == 5
        plot(squeeze(perturb_act(site,5,:)),'k')
    else
        plot(mean(squeeze(perturb_act(site,5:end_points(site),:)),1),'k')
    end
    title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end


% average tp - running onsets
end_points = [5, 6, 9, 9, 9, 10];
figure;
for site = 1:size(proj_meta,2) 
subplot(ceil(size(proj_meta,2)/2),2,site);
hold on
plot(mean(squeeze(run_act(site,1:2,:)))','b')
plot(mean(squeeze(run_act(site,3:4,:)))','r')
if end_points(site) == 5
    plot(squeeze(run_act(site,5,:)),'k')
else
plot(mean(squeeze(run_act(site,5:end_points(site),:)),1),'k')
end
title([regexprep(proj_meta(site).animal,'\_','\\_') ' - ' num2str(proj_meta(site).ExpGroup(1))]);
end

% average tp - mismatch - across all animals
end_points = [5, 6, 9, 9, 9, 10];
dim = size(perturb_act);
temp6=squeeze(perturb_act(1,5:5,:))';
for ind=2:dim(1)
    temp6 = [temp6; squeeze(perturb_act(ind,5:end_points(ind),:))]; 
end

aa=figure;hold on
ac=(mean(reshape(perturb_act(:,1:2,:),2*dim(1),dim(3)))-1)*100;
ac_sem=SEM(reshape(perturb_act(:,1:2,:),2*dim(1),dim(3)))*100;
ac2=(mean(reshape(perturb_act(:,3:4,:),2*dim(1),dim(3)))-1)*100;
ac_sem2=SEM(reshape(perturb_act(:,3:4,:),2*dim(1),dim(3)))*100;
ac3=(mean(temp6)-1)*100;
ac_sem3=SEM(temp6)*100;
area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac + ac_sem fliplr(ac-ac_sem)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac2 + ac_sem2 fliplr(ac2-ac_sem2)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac3 + ac_sem3 fliplr(ac3-ac_sem3)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');

plot(-l_win/10:1/10:r_win/10,ac,'Color','b','Linewidth',3)
plot(-l_win/10:1/10:r_win/10,ac2,'Color','r','Linewidth',3)
plot(-l_win/10:1/10:r_win/10,ac3,'Color','k','Linewidth',3)

xlim([-1.1 5.2])
ylim([-2 7])
xlabel('time [s]')
ylabel('\DeltaF/F [%]')

%%% same stuff all post lesion pooled
% average tp - mismatch - across all animals
end_points = [5, 6, 9, 9, 9, 10];
dim = size(perturb_act);
temp6=[]
for ind=1:dim(1)
    temp6 = [temp6; squeeze(perturb_act(ind,3:end_points(ind),:))]; 
end

aa=figure;hold on
ac=(mean(reshape(perturb_act(:,1:2,:),2*dim(1),dim(3)))-1)*100;
ac_sem=SEM(reshape(perturb_act(:,1:2,:),2*dim(1),dim(3)))*100;
% ac2=(mean(reshape(perturb_act(:,3:4,:),2*dim(1),dim(3)))-1)*100;
% ac_sem2=SEM(reshape(perturb_act(:,3:4,:),2*dim(1),dim(3)))*100;
ac3=(mean(temp6)-1)*100;
ac_sem3=SEM(temp6)*100;
area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac + ac_sem fliplr(ac-ac_sem)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
% area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac2 + ac_sem2 fliplr(ac2-ac_sem2)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');
area([-l_win/10:1/10:r_win/10 r_win/10:-1/10:-l_win/10],[ac3 + ac_sem3 fliplr(ac3-ac_sem3)],'Facecolor',[0.9 0.9 0.9],'LineStyle','none');

plot(-l_win/10:1/10:r_win/10,ac,'Color','b','Linewidth',3)
% plot(-l_win/10:1/10:r_win/10,ac2,'Color','r','Linewidth',3)
plot(-l_win/10:1/10:r_win/10,ac3,'Color','k','Linewidth',3)

% xlim([-1.1 5.2])
% ylim([-2 7])
xlabel('time [s]')
ylabel('\DeltaF/F [%]')
figure,hold on
for ind=1:4
%     plot(mean(squeeze(perturb_act(:,ind,10:end)))')
ar(ind) = max(mean(squeeze(perturb_act(:,ind,10:30)))')
% % pause
end
for ind=5:9
    plot(mean(squeeze(perturb_act(3:end,ind,1:end)))')
ar(ind-2) = max(mean(squeeze(perturb_act(3:end,ind,10:30)))');
% figure,plot(squeeze(perturb_act(3:end,ind,1:end))')
% pause
end

ar(1)=max(mean(reshape(perturb_act(:,1:2,:),2*dim(1),dim(3))));
ar(2)=max(mean(reshape(perturb_act(:,3:4,:),2*dim(1),dim(3))));
